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ABSTRACT 

The formation of pillars of dense gas at the boundaries of H II regions is investigated with 
hydrodynamical numerical simulations including ionising radiation from a point source. We 
show that shadowing of ionising radiation by an inhomogeneous density field is capable of 
forming so-called elephant trunks (pillars of dense gas as in e.g. M16) without the assistance 
of self-gravity, or of ionisation front and cooling instabilities. A large simulation of a density 
field containing randomly generated clumps of gas is shown to naturally generate elephant 
trunks with certain clump configurations. These configurations are simulated in isolation and 
analysed in detail to show the formation mechanism and determine possible observational 
signatures. Pillars formed by the shadowing mechanism are shown to have rather different 
velocity profiles depending on the initial gas configuration, but asymmetries mean that the 
profiles also vary significantly with perspective, limiting their ability to discriminate between 
formation scenarios. Neutral and molecular gas cooling are shown to have a strong effect on 
these results. 
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1 INTRODUCTION 

Among the most striking features of H II regions are the long 
columns of neutral gas pointing towards the central star, known 
variously as elephant trunks (ETs), pillars and fingers. The most 
famous examples are the 'pilla rs of creation' in th e Eagle Nebula 
(Ml 6), observed with HST bv lHester et all dl996h . The HST im- 
ages of NGC 3372 (the Carina Nebula ) clearly show a num ber of 
ETs, as do HST images of NGC 3603 terandner et ai]|2000l). The 



Elephant Trunk Nebula in IC 1396, obser ved by e.g. Reach et al 



2004) is another typical example. ICarlqvist. Gahm. & Kristen 



20031) show ETs in the Rosette Nebula (NGC 2237-2246), IC 
1805, and N81. 

The Ml 6 pillars have been well studied observationally. They 
are towards the more massive end of the scale of observed ETs but 
are not extreme, so we ha ve used their proper ties as a benchmark 
for our simulation results. Hest er et"al1 d 19961) studied the surface 
of the M16 pillars in detail with HST, finding that the interface be- 
tween the pillars and the H II region is well explained by a thin ion- 
isation front (I-front) with a strong photo-evaporation flow into the 
lower density ionised gas. Due to the large opacity of the molecu- 
lar gas they could say little about the interior of the pillars. Molec- 
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ular CO maps were obtained by IPounoj dl998l) . who showed that 
the pillars are mostly molecular, that the internal and external pres- 
sures are comparable, and that their photo-evaporation times are 
about 2 x 10 7 years. He also found strong velocity gradients in the 
molecular emission, but it is not cl ear if this indicates a coherent 
flow generated by pillar formation. Whi te et al. I dl999l) . using in- 
frared, millimetre and radio observations, found the pillars to have 
cold cores (~ 20 K), with surrounding warm (~ 60 K) gas, and an 
outer 'hot' shell (at ~ 250 — 320 K). The clumps at the heads of the 
pillars appeared to have typical densities 2 x 10 s cm -3 , while the 
trunks exhibited lower mean densities of 3 — 5 x 10 4 cm" 3 . These 
authors also found that the internal pressure of the heads of the pil- 
lars is a factor of ~ 2 lower than the adjacent ionised gas pressure. 
They interpreted this as evidence that the heads o f the pillars ar e 
in the early stages of radiatively driven implosion ( lBertoldil ll989), 
an d may therefore be < 100 kyr old. Subsequent dy namical models 
by Williams, Ward-Thompson, & Whitworthl d200ll) suggest alter- 
nate scenarios in which the pillars could be 200 — 400 kyr old. 

Despite detailed observations, and a long history of theoretical 
models, there is as yet no clear consensus as to how such features 
as ETs arise, their lifetimes, or their end states. It has been sug- 
gested that instabilities in I-fronts can generate similar structures or, 
alternatively, that pre-existing inhomogeneities in the interstellar 
medium (ISM) will create shadowed regions behind dense clumps 
where gas can accumulate. The question of which of these s cenar- 
ios ac tually forms ETs was discussed at least as far back as Kahn 
dl958l) . Given that the ISM is clumpy, and that instabilities are 
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present under certain conditions, it is likely that both processes con- 
tribute to some degree. Thus the difficulty from a theoretical per- 
spective is that there are a number of potentially viable mechanisms 
to form pillars, and it is difficult to observationally di stinguish 
betwe en them. The instability model was proposed by iFriemarj 
d 1954) for Rayleigh-Taylor instabilities, and subsequently devel- 
oped for I-front instabilities. These instabilities have been studied 
by many authors but are not the s ubject of this pa p er, so we re- 
fer the reader to recent work by e. g. lWilliams! J20 02); Mizut aet al.l 
(2006): lWhalen & Normar] (2008) and references therein. We con- 
centrate here on modelling the effects of the sha dowing of dense 
clump s on the dynamics of an expanding I-front. iLim & Mellemal 
(2003) found that a partially shadowed clump in an ionising ra- 
diation field was pushed further into the shadowed region, and a 
denser and longer-lived tail arose compared to that of a single iso- 
lated clump. T his was in the context of Earth-mass clumps in plan- 
etary nebulae dO'Dell & Handronll 19961) which have lifetimes of a 
few thousand years. For more massive clumps, the longer lifetimes 
may allow much denser tails to develop in the shadowed regions, 
possibly to the extent that they would be observed as pillars. This 
was also suggested bv lPound et al.l d2007l) . 

General models of the evolution of a photo-evaporating 
clu mp and the shadowed regio n behind it were devel oped 
by iBertoldi & McKed jl990t) and iLefloch & Lazaref j dl994l) . al- 
though these were m ostly concerned with the dense clump itself. 
Willi ams et all feOOlb used computational models to study forma- 
tion scenarios of the Ml 6 pillars, setting up axi-symmetric simu- 
lations with various initial density fields exposed to planar ionis- 
ing radiation. They found that with their modelling assumptions, 
multiple initial scenarios were capable of producing dense struc- 
tures resembling the ETs seen in M16. These were found to be 
long-lived quasi-equilibrium structures, raising the possibility that 
the pillars could be quite old (> 300 kyr). Their work highlights 
the difficulty in interpreting the observations - a number of differ- 
ent initial con ditions and physica l processes could produce pillars. 
More recently iMiao et al.l yOOfj) used an SPH code with radiative 
transfer to model the photo-ionisation of a dense cloud, their re- 
sults supported the idea that the head of the pillar is in the early 
stages of radiation-driven implosion. Their modelling was of the 
head more than the trunk, however, and it remains unclear whether 
a pillar could fo rm behind the im ploding head on a time-scale as 
short as 100 kvr. lPound et alj J2007h modelled the photo-ionisation 
of a single dense clump (M — 30 Mq) by an O star, specifically 
looking at the shadowed region. They found the shadow can pro- 
duce a long neutral tail of the dimensions of the M16 pillars, but 
without enough material in the tail. They note that adding in mul- 
tiple clumps of different sizes should increase t he amount of tail 
material, as was found bv lLim & Mellemal (2003), a suggestion we 
explore in detail in this paper. 

The first global simulation of the expansion of an H II re- 
gion in to a turbulent density field was presented bv lMellema et al.l 
d2006al) . They found that features such as ETs developed quite 
naturally due to the uneven I-front expansion velocity, but due 
to the nature of the simulati on individual pillars are poorly re- 
solved. iMac Low et al. I d2007l) also modelled global expansion of 
an H II region, but their simulations were more of the early I- 
front e xpansion than of the late r dynamical evolution. Very re- 
cently, iGritschneder et all d2009l) modelled part of an expanding 
H II region with planar radiation impinging on a turbulent density 
field. They also found that pillar-like features developed naturally 
in their simulations after about 250 kyr, but again t he res olution in 
individual pillars is poor. lLora. Raga. & Esquivell j2009h also fol- 



low a similar approach to form pillar-like structures. They study 
the angular momentum of dense clumps which form in their simu- 
lation and find preferential alignment perpendicular to the direction 
of the radiation field. In wo rk that is in some ways similar to ours, 
but on much smaller scales, iRaga et all (2009) showed how photo- 
evaporation flows from a large reservoir of dense gas can flow into 
shadowed regions, recombine, cool, and begin to build up dense 

pillar-like structures. 

With some excep tions dWilliams et al.l 1200 ll ; IPound et al.l 
l2007l ; lRaga et alj |2009), these works were not primarily focussed 
on how ETs form, and as a result it is difficult to tell what phys- 
ical processes are most relevant. The aim of this work is to focus 
on the shadowing mechanism to see how effectively it can produce 
ETs in an idealised environment. We have developed a radiation- 
magnetohydrodynamics code with which to study this problem. We 
will describe our code and algorithms in section|2] In section[3]we 
describe the initial conditions and show results from 3D simula- 
tions of the photo-ionisation of a density field with randomly dis- 
tributed dense clumps. These models show a number of structures 
resembling ETs which form dynamically due to shadowing in the 
inhomogeneous medium. In section [4] we simulate certain clump 
configurations in isolation to demonstrate two different ways pil- 
lars could form. The first model has the clumps oriented almost 
like a pillar in the initial conditions, whereas in the second model 
clumps are swept past each other into a pillar-like structure. We 
find that neutral gas cooling has a strong effect on our results and 
in section [5] we repeat these two models using an alternate ther- 
mal model with moderately strong neutral gas cooling. The details 
and evolutionary time-scales change considerably, but the forma- 
tion mechanisms remain unchanged. In section [6] we discuss the 
context and significance of our results and in section [7] we present 
our conclusions. 



2 NUMERICAL METHODS AND ALGORITHMS 
2.1 Fluid Dynamics 

We have written a modular, finite volume, fluid dynamics code to 
run these simulations. The code uses a uniform grid in 1, 2, or 3 
spatial dimensions with cubic cells. Scalar and vector fields are 
both cell-centred. The integr ator for the fluid equations is based o n 
the algorithm described by iFalle. Komissarov. & Joarderl dl998l) . 
which is dimensionally unsplit and second order accurate in time 
and space. We have separate Riemann solvers for the Euler and 
Ideal Magnetohydrodynamics (MHD) equat ions. We also add some 
artificial viscosity in a similar manner to Falle et al.l d 19981) us- 
ing a coefficient of rj v — 0.05 — 0.1. This is required to ensure 
shocks travel at the correct speed in all directions. It fixes for ex- 
ample the "carbuncle proble m" in the Double Mach Reflection test 
dWoodward & Colellal 1981 used the Lapidus viscosity pre scription 
for this) and mitigates the Quirk instability jOuirkl fl994l) for sta- 
tionary grid-aligned shocks. 



2.2 Ray-tracing and Microphysics 

Our ray-trac ing and microphysics r outines are based largely on the 
methods in Lim & Mellema (2003) and on the C 2 -ray method de- 
veloped by Mell ema et al.l ( l2006bl) . We use operator splitting to first 
update the dynamics by a full timestep, and then run the micro- 
physics update over the full timestep. In this work we only consider 
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explicitly the ionisation of atomic hydrogen. We first describe the 
ray-tracing algorithm and then the microp hysics calculation . 

The Short Characteristics tracer (e.g. lRaga et al.ll999T) is used 
to trace out rays from a source in a causal manner, calculating the 
optical depth to a cell by interpolating between (previously cal- 
culated) optical depths to neighbouring cells closer to the source. 
Given that we are ignoring diffuse radiation (the On-the-Spot ap- 
proximation) the diffusion in the ray-tracer is not significant, and 
is minim ized using the weighting scheme given bv lMellema et al.l 
d2006bl) . 

When the photo-ionisation time is short compared to other 
time-scales (cooling, recombination, and collisional ionisation 
times) the microphysics equations become difficult to solve explic- 
itly so we adopt a dual approach. In cases of weak photo-ionisation, 
we use an explicit 5th order Runge-Kutta technique with adap- 
tive step-size to a given relative accuracy jPress et alJl 19921) . For 
strong photo-ionisation we integrate explicitly until the Hydrogen 
ion fraction, x, satisfies x 0.95, and then analytically integrate 
th e equations assuming a constant electron density (as described 
in llVlellema et al.l l2006bh . with bisection substepping to conver- 
gence (typically 2 — 4 substeps). For both of these methods we 
use a relative error tolerance of 0.001. 

This algorithm also calculates the time-averaged optical depth 
through the c ell At, which i s then used by subsequent cells in 
the ray-tracer. iMellema et ail {2006b) use a simple time average 
of At, however we use a time average of exp(— At) since this 
gives a time average of the fraction of photons passing through the 
cell. This can be easily seen in the (extreme) case of an optically 
thick cell which is photo-ionised "rapidly" half way through a unit 
timestep, so that 



Ar(t) = 



100 t < 0.5 
0.5 < t < 1 



(1) 



The mean optical depth over the timestep is 50, but clearly 
half of the incident photons will pass through the cell, and 
Jq 1 exp(— At)(U = 0.5 gives the desired result. We do this inte- 
gration at the same time as the microphysics variables, to the same 
accuracy criterion. 

We use monochromatic radiation with a hydrogen photo- 
ionisation cross-section of 6.3 x 10 -18 cm 2 and an energy gain of 
5.0 eV per photo-ionisation. Col lisional ionisati on rates are calcu- 
lated with fitting functions from lVoronovl dl997l) . and ra diative re- 
combi nation (Case B) rates using the tables calculated bv lHummeil 
dl994l) . The difference between planar radiation and radiation from 
a point source can be quite significant if the size of the compu- 
tational domain is comparable to the distance to the source. The 
rocket effect is weaker further from a point source due to the in- 
verse square law, which may extend the lifetime of any structures 
that form. This effect can, however, reduce the length of such struc- 
tures since the intensity of the radiation is higher at their heads. In 
the case of Ml 6, the heads of the pillars are about 2pc from the 
brightest star, and they are about 1 pc long, so the flux dilution is 
more than a factor of 2 along their length. We therefore use a point 
source in this work. 



2.3 Gas Cooling 

We use two cooling models in this work, denoted CI and C2, which 
differ in that C2 has significant cooling in the neutral gas. Our first 
model, CI, contains four elements: 

(i) Radiative losses due to recombining Hydrogen, calculated 



from the non-equilibrium ion fraction and tem perature in each cell 
according to rates tabulated in lHummej ( fl99i) . 

(ii) Collisional ionisation of Hydrogen: this is relatively unim- 
portant because the rates are typically very low, but we subtract the 
ionisation energy from the gas for each collisional ionisation. 

(iii) Cooling due to heavy elements at high temperatures, using 
th e collisional ionisation equ ilibrium (CIE) cooling curve tabulated 
in lSutherlan d & Dopita| { 19931) and shown in their fig. 18. This pro- 
vides strong cooling in ionised gas with temperatures significantly 
larger than 10 000 K. Note that in CIE at 10 000 K, Hydrogen, Ni- 
trogen and Oxygen are neutral so this fitting function does not dou- 
ble count the other terms in our cooling function, at least for the 
gas temperatures encountered in our simulations. 

(iv) A linear fit to collisionall y excited emi ssion from photo- 
ionised Oxygen and Nitrogen dOsterbrockl l 1989). 

The last term is the most important for this work, since these ionic 
species are the dominant coolants in H II regions and set the equi- 
librium temperature in ionised gas of T eq ~ 8000 K. In experi- 
ments with different cooling functions for ionised gas, we found 
that the most important factor for the dynamical evolution of our 
models was the equilibrium temperature. If the normalisation of 
the cooling function is kept fixed at 8000 K, its slope has little ef- 
fect on the resulting dynamics so long as the slope is positive. If 
we had strong shocks in the ionised gas this aspect of the cooling 
function would have more influence, but the photo-ionised gas in 
our simulations has a very narrow temperature range. 

In this prescription the neutral atomic gas has no efficient cool- 
ing avenue, and shocked neutral gas is typically at 100 — 10 000 K. 
This is undoubtedly a limitation in our modelling, but we do not 
yet model the formation of molecules, or the formation/destruction 
of dust, which are the primary neutral gas coolants in star forming 
regions. To assess the effects of significant neutral gas cooling we 
also use an alternate cooling function, C2, consisting of the previ- 
ous components in CI plus additional exponential cooling (New- 
ton's Law) in neutral gas with a rate given by 



T = 



(1 



lO^yrs 



(r^ - T) 



(2) 



where T is gas temperature, Too is the temperature to which this 
cooling law relaxes at late times, x is the ionisation fraction of 
the gas, and TV is a parameter specifying the chosen cooling time- 
scale, t c = 10 Jv yrs/(l — x) 2 . The scaling with (1 — x) 2 en- 
sures only mostly neutral gas is affected. We set Too = 100 K 
and TV = 4 for the alternate models run in this paper. This is 
not an extreme model either in terms of the equilibrium temper- 
ature or the cooling tim e, having less coolin g in dense gas than 
the model presented in iHennev et al. I d2009h . It is a very sim- 
ple prescription, with an effect which is inte rmediate between CI 
and a two-temperature isothermal model ( e.g. lWilliams et alj200ll : 
iGritschneder et alj2009l : lLora et alj|2009l) . 



2.4 Code Tests 

We have extensively tested the fluid dynamics, microphysics, and 
raytracing components of the codeQ ] For hydrod ynamics we used 
a range of shock- tube tests in ID dTorol il999), and then in 2D 
at various angles to the grid axes, with the code reproducing 
the correct solutions. We have run the double mach reflection 



1 A brief description of our code and results from test problems can be 
found at http : / /homepages . dias . ie/ " jmackey / jmac/ 
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Figure 1. Tests of the radiative transfer algorithm in 2D (left) and 3D (right). The top two plots compare the number of ions to the number of photons emitted 
for a source in a uniform medium with dynamics and recombinations switched off. The source is at the centre of the domain, and the I-front remains on the 
domain for the duration of the simulations. This is a scale-free problem which we have run with parameters such that the cell optical depth is St = 10 for the 
runs with 101 2 and 101 3 cells (St = 30.6 for lower resolution and 3.93 for higher resolution). We always lose some photons due to the interpolation, but this 
decreases dramatically with resolution. The lower plots show the evolution of the I-front radius over time for runs with recombinations switched on. Again it 
is a uniform medium with the same optical depths per cell. The time is shown in units of the recombination time t rcc , and the radius in units of the Stromgren 
radius (or its 2D analogue). When the timestep St = t Iac the I-front propagates too slowly, but for St = O.ltrec it has the correct speed. 



test dWoodward & Colellal |l984) and obtained good agreement 
with t he original work and with e.g. the ATHENA code dStone et all 
fig. 16). Since we explicitly add in numerical viscosity, the 
diffusion is slightly stronger than in the ATHENA code. We have 
also done tests of implosions and blast wav es in 2D and 3D, fi nding 
results consistent with previous work (e.g. IStone et al . 2008), and 
recovering the Sedov-Taylor solution for the adiabatic blast wave in 
3D. We have a lso tested the deve lopment of the Kelvin-Helmholz 
instability (e.g. lAgertz et all2007l) . finding very satisfactory agree- 
ment with other work. 

iRavmoncj |j979) calculated zero-dimensional shock models 
by following a parcel of gas through a shock front. We have tested 
our code against his 'Model E', a 100 km s -1 steady shock with 
a 1/iG transverse magnetic field, setting it up in our code as a 
ID problem with 100kms -1 gas hitting a dense cold layer and 
allowing the system to relax to an equilibrium state. Our test re- 
produced well the ion-fractions for Hydrogen and Helium, the 
gas temperature and the density as a function of position. We 
have also run models at higher shock velocity and find that for 
v > 130 km s _1 the shocks become overstable, in agreement with 
previous work (e.g. llnnes. Giddings, & Falle 19§3). 

Photo-ionisation was tested in conjunction with raytracing us- 



ing similar tests to those in iMellema et al. (2006b), where the dy- 
namics is switched off. We started with ID rays from a source 
at infinity, without dynamics or recombinations. For a grid with 
1 000 cells, we computed models with cell optical depths At = 
{0.1, 1, 10, 100}, and where the total number of timesteps varied 
from Uim/St = { 10 1 , 10 2 , 10 3 , 10 4 } . The error in I-front position 
compared to the analytic value was found to converge rapidly to 
less than one cell width with increasing time resolution. For mod- 
els with recombinations turned on, errors were no more than than 
one cell width for all runs with > 10 timesteps per recombination 
time, except for low density models where the I-front is resolved. 

In 2D and 3D, we computed the expansion of circular and 
spherical I-fronts from a point source into a static medium, with 
and without recombinations. Without recombinations, the models 
provide a test of photon conservation (by comparing the number of 
ions to photons emitted as a function of time). We plot the pho- 
ton conservation in the top two panels of Fig. [T] for 2D on the 
left and 3D on the right. These figures show the relative sizes of 
ray-tracing and time-integration errors as a function of resolution 
and dimensionality. For the very low resolution runs (33 2 and 33 3 
cells) we lose between 1 and 10 per cent of photons due to inter- 
polation errors in the ray-tracing when the ionised region is < 10 
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cells across. With increased spatial resolution the errors decrease 
strongly whereas increased time resolution doesn't help the 33 cell 
runs significantly. There is a dramatic improvement in accuracy 
with time resolution for the 101 and 257 cell runs. These results 
show that the errors are interpolation dominated when the number 
of cells is much smaller than the number of timesteps and time- 
integration dominated in t he opposite limit. Usin g the weighting 
scheme recommended by Melle ma et ail (2006b) we find that I- 
fronts are circular to within a cell width over a wide range of den- 
sities, luminosities, and spatial and temporal resolutions. 

The lower panels of Fig. Q] show the position of the I-front 
as a function of time for simulations with recombinations in- 
cluded, modelling the idealised Stromgren Sphere analysis and its 
2D analogue. The mean I-front position is always within 1 — 2 
per cent of the analytic value except at very early times when it 
has only crossed a few cells, or when the timesteps are of or- 
der the recombination time, t TCC = (anii) -1 , where a is the 
(case B) recombination coefficient set to a constant for this test 
(a = 2.59 x 10 -13 cm 3 s _1 ), and tih is the Hydrogen number 
density. This is an expected limitation of the C 2 -ray method since 
it uses time-averages of the photon flux through each cell (see 
Melle ma et alj|20 06b). For the tests where St — t rcc we underesti- 
mate the I-front velocity while it expands to the Stromgren radius. 
The error is slightly larger at higher spatial resolution because we 
have to do the same inaccurate time-average across more cells and 
the error is always on the side of losing photons. For sufficient time 
resolution, however, the I-front propagates at the correct speed, and 
it is worth noting that the photo-ionisation time for a cell is much 
shorter than the recombination time while the I-front is expanding 
rapidly. We do not need to resolve this timescale to get accurate 
results. This is the major strength of the C 2 -ray algorithm. 

We have also used our code to reproduce the simulations 
of iLim & Mellemi] J20031) . where dense clumps in different con- 
figurations were photo-evaporated by planar ionising radiation. We 
find our code gives largely consistent result s despite significant dif - 
ferences in the numerical methods, e.g. the lLim & Mellemal J2003h 
simulations did not use an explicitly photon-conserving algorithm. 



3 RANDOM CLUMPS SIMULATIONS AND RESULTS 



Previous studies (e.g. Mellema et al. 2006a; Ma c Low et al. I l2007l ; 
iGritschneder et all2009h have used a turbulence model to generate 
a density field into which the ionising radiation propagates. In this 
work we wish to isolate the shadowing effect of the radiation from 
other rather uncertain physics such as the type of turbulent motions 
we inj ect into the gas. In thi s respect our work is more similar to 
that of Williams et "all fcOOlh . We add dense clumps to a uniform 
density field in the following way: 

(i) We start by setting the mean background density, rib, ranging 
from 10 to 10 3 cm -3 in different simulations. 

(ii) We choose a total mass to put into clumps by calculating the 
mass associated with a smoothed mean density in clumpy material, 
rid, in a subset of the full simulation domain. The subset is chosen 
to keep clumps away from simulation boundaries, particularly the 
boundary nearest the radiation source. 

(iii) We set the maximum and minimum mass and size of the 
clumps and draw clumps randomly within these limits until all the 
mass in clumpy material has been used up. The size of a clump 
in each direction is chosen separately, so typically we get triaxial 
clumps. 



(iv) The clump's mass and size determine its peak overdensity 
for a given radial profile. In this work we use Gaussian profiles. 

(v) The clump's position is then set at random on the domain, 
within the allowed region, and it is rotated by a random angle in all 
three directions. Clumps are added to the background density field 
one by one, in a linear superposition when they overlap. 

(vi) The initial conditions are static everywhere, and we assign 
a constant pressure throughout. 

The method of random sampling can strongly influence the 
distribution obtained. Clump positions, radii, orientations, and 
masses are selected randomly within certain minimum and maxi- 
mum values. In all cases, we select a random number on a linear 
scale between the two limits. 



3.1 2D Simulations 

The above procedure can also be performed in 2D with slab sym- 
metry, enabling us to do a study with a wider range of parame- 
ters than is possible in 3D. We performed 184 2D simulations on 
a 1.5 x 1.5 pc domain with a radiation source 2pc off the do- 
main in one direction. The background densities used were rib = 
{10, 100, 1000} cm -3 ; the ionising photon fluxes entering the do- 
main were F-, = {10 10 , 10 11 , 10 12 , 10 13 } cm* 2 s" 1 ; the num- 
ber of clumps were N c i — {10,30, 100}, being either compact 
or extended, round or triaxial, with scale radii 0.0225 pc ^ R s ^ 
0.15 pc. Compact round clumps had radius 0.05 pc, extended round 
clumps were 0.12 pc, compact triaxial clumps were [0.0225 — 
0.06] pc, and extended triaxial clumps were [0.06 — 0.15] pc. The 
total mass in clumps was chosen by the above method, setting a 
smoothed mean density in clumps of n c ; = {10 2 , 10 3 , 10 4 } cm -3 
and calculating the associated mass. For the 2D models we set the 
maximum and minimum clump masses to be very close to each 
other so that all clumps would have similar masses. 

In addition to the 184 models with the parameters mentioned 
above, we also selected particular models and varied the random 
seeds, also using larger domains to ensure boundary effects were 
not important. For models with 10 clumps, the specific clump con- 
figuration was important, but for 30 and 100 clumps, the random 
seed had little effect on whether or not the model produced ETs. 
The boundary effects were found to be minimal. We use outflow, or 
zero-gradient, boundary conditions for all of the simulations in this 
paper. This is the logical choice for the lateral boundaries and for 
the boundary furthest from the source. For the boundary nearest the 
source a reflection boundary could be imposed to model a pressure 
confined H II region, but we have chosen to model an H II region 
which has broken out of its parent molecular cloud and is no longer 
pressure confined in at least some directions. 

Care must be taken with the transition from R-type to D- 
type I-fronts. If this occurs too close to the boundary nearest the 
source spurious inflows can be set up by rarefaction waves prop- 
agating backwards from the stalling I-front. This can introduce an 
unphysical aspect to the simulation results, therefore we chose all 
of our simulation parameters such that the I-front in the background 
medium remained R-type well into the domain. This ensured that 
photo-evaporation flows dictated the flow at the boundary rather 
than vice versa. We have verified that the R-type to D-type transi- 
tion occurs as expected with a shock driven ahead of the I-front into 
the neutral gas and a strong photo-evaporation flow established in 
the opposite direction. 

Of the 184 2D models, approximately 30 produced pillar-like 
structures, with another ~ 20 producing short features more like 
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No. 


n c i 




M cl 


N cl 


Size 


Mi(min/max) 






Mj 


Flux 


Results 


1 


10 3 




72.2 


115 


C 


0.11 


- 1.1 


2.7 x 


10 s 


12.4 


10 10 


Nothing like ETs; corrugated I-front. 


2 


10 3 




72.2 


115 


c 


0.11 


- 1.1 


2.7 x 


10 s 


12.4 


10 11 


Nothing like ETs; corrugated I-front. 


3 


10 3 




72.2 


103 


E 


0.36 


- 1.1 


2.0 x 


10 4 


45.3 


10 11 


Too diffuse to form dense structures. 


4 


10 3 




72.2 


40 


C 


0.36 


- 2.9 


9.1 x 


10 s 


6.7 


10 10 


Low mass and low density ETs form. 


5 


10 3 




72.2 


40 


C 


0.36 


- 2.9 


9.1 x 


10 s 


6.7 


10 11 


Short lived Cometary structures. 


6 


10 3 




72.2 


40 


E 


1.3 - 


- 2.6 


4.1 x 


10 4 


31.7 


10 11 


No dense ET-like structures. 


7 


10 3 




72.2 


54 


C 


0.22 


- 2.2 


4.8 x 


10 s 


9.3 


10 10 


Nothing dense enough or trunk-like. 


8 


10 3 




72.2 


54 


C 


0.22 


- 2.2 


4.8 x 


10 s 


9.3 


10 11 


Nothing dense enough or trunk-like. 


9 


10 3 




72.2 


85 


C 


0.14 


- 1.4 


4.0 x 


10 s 


10.1 


10 10 


Small ET-like structures. 


10 


10 3 




72.2 


85 


C 


0.14 


- 1.4 


4.0 x 


10 s 


10.1 


10 11 


Small cometary structures, low density. 


1 1 


10 4 




666 


118 


C 


1.0 - 


- 10 


3.5 x 


10 e 


3.4 


10 11 


Some ETs among large 'mountains'. 


12 


10 4 




666 


45 


C 


3.3 - 


- 27 


7.6 x 


10" 


2.3 


10 11 


Multiple dense long-lived ETs. 


13 


10 4 




666 


40 


E 


12 - 


24 


5.7 x 


10 s 


8.5 


10 11 


Large 'mountains' with a few dense ETs. 


14 


10 4 




666 


59 


C 


2.0 - 


- 20 


5.1 x 


10 e 


2.8 


10 11 


Many dense ETs. 


15 


10 4 




666 


83 


C 


1.3 - 


- 13 


3.9 x 


10 6 


3.3 


10 11 


Many dense ETs. 


16 


4 x 10 3 


266 


40 


E 


2.4- 


- 10.4 


1.8 x 


10 s 


15.1 


10 11 


Formed long-lived ETs, see Fig. [2] 


17 


n b = 


200 


84.9 


3 


0.1 pc 


28.3 




1.0 x 


10 s 


20.1 


2.9 x 10 11 


Almost co-linear clumps, see Fig. [3] 


18 


n b = 


200 


83.7 


3 


0.1 pc 


27.7 


- 28.3 


1.0 x 


10 s 


20.1 


2.9 x 10 11 


Triangular configuration, see Fig. [5] 



Table 1. Simulation parameters for the test 3D simulations (1 — 15) and the one high resolution model (16). All models have a background density n b = 
100 cm -3 and triaxial clumps as defined in the text, except for models 17 and 18 which have m, = 200 cm -3 and spherical clumps (see section|4j. Clump 
density, n c ; refers to the mean density of gas put into clumps ( cm -3 ) and the next column M c ; gives the total mass in clumps (all masses are in solar masses, 
Mq). N c i lists the number of clumps in each model. For size, 'C refers to compact and 'E' to extended, as described in the text (section l3~TV Mi(min/max) 
gives the range of clump masses for the model. The Jeans mass (defined as Mj = 2O3M0 T 15 n — 0,5 for adiabatic gas) is calculated for the densest point 
in the initial conditions (n max ) and for a temperature of 10 K. We do not use the actual clump temperatures because the initial constant pressure state has 
very low temperatures at the densest points (T ~ 1 K in some cases); this gives a meaningless Jeans mass since if we included gravity we wouldn't need to 
impose such an artificially low temperature. The Flux is the monochromatic ionising photon flux (cm" 2 s" 1 ) entering the domain. All models use clumps 
with random positions and properties (within ranges) except models 17 and 18 which are set up to model specific clump configurations in isolation. 



heaps than pillars, and another ~ 10 producing very short-lived 
pillar-like features lasting < 50kyr. The most successful groups 
of models had background and clump densities of [n b ,n c i] — 
[10 2 ,10 3 ]cm" 3 , [10 2 ,10 4 ]cm" 3 , and [10 3 , 10 4 ] cm" 3 , with 10 
or 30 clumps. Those models which successfully produced pillars 
had a combination of large low-density regions where the I-front 
could propagate far from the source, and large dense clumps which 
were capable of stalling the I-front for sufficient time to produce a 
pillar in the shadowed region. The trunks of the pillars were formed 
from clump material more than ambient gas, simply because the 
ambient gas is not dense enough. If the background density was 
high enough to contribute significantly to a pillar, then it was also 
high enough to stall the I-front and a pillar could not form. The 
clump material which made up a pillar always came from several 
clumps; a lone clump in these models never gathered enough ambi- 
ent gas into the tail to produce a pillar, nor did it contribute enough 
of its own gas to the shadowed region to form the trunk of a pillar. 
Typically, pillars formed via the interaction of a number of clumps 
shadowing each other and accelerating past/through each other pro- 
duced a dense tail which lasted anything from 50—200 kyr, depend- 
ing on the inertia of the leading clump. 

3.2 3D Test Simulations 

There are significant differences between the 2D and 3D cases: the 
stellar flux drops off more rapidly with distance in 3D; shocks in 
the shadowed regions are converging or diverging whereas in 2D 
they are planar; there is also an extra dimension for gas to flow 
past a clump. These favour the formation of ET-like features in 
2D, however our 2D results can be used to constrain our choice 
of parameters for more computationally expensive 3D simulations. 
Some of these differences can be modelled by comparing 2D slab- 



symmetric and axi-symmetric models. For our models, however, 
it is the asymmetries which generate the ETs and these cannot be 
modelled in axi-symmetry. 

Initially we computed fifteen 3D models with 168 x 128 2 
grid cells, since the 2D models had shown that this was sufficient 
to resolve the larger features in the simulations. The physical 
domain was 2.625 x 2 x 2pc 3 . The parameters we varied in 
each model are listed in Table [T] models 1 — 15. All models start 
with a constant initial pressure of Pg/ks = 2 x 10 5 cm -3 K 
where ks = 1.38 x 10 _16 ergK _1 is the Boltzmann constant; 
this pressure is higher than the typical ISM value and is a com- 
promise between not having the temperature too high in diffuse 
gas and too low in dense gas. The results are not sensitive to 
this initial pressure because the ionised gas pressure is 100 x 
greater and it is the pressure difference which drives the dynamical 
evolution of our models. All 3D models had n b = 10 2 cm -3 ; 
of these four had {n c! ,F 7 } = {10 3 cm" 3 , 10 10 cm' 2 s" 1 }, 
six used {10 3 cm" 3 , 10 11 cm" 2 s" 1 }, and five used 
{10 4 cm" 3 , 10 11 cm" 2 s _1 }. Of these, the models with 
n c i = 10 3 cm -3 had too little mass in clumps and none of 
the models produced anything as dense or massive as a pillar. 
The models with n c i — 10 4 cm" 3 , by contrast, did produce 
massive and long-lived pillars, simply because the initial clumps 
were so massive and dense. The models with fewer clumps were 
more successful. With too many clumps there were no large 
low-density regions, resulting in structures more like mountains 
than long columns. These dense models were rather unrealistic 
however, with initial densities in clumps of «h > 10 6 cm -3 , 
very gravitationally unstable, as can be seen by comparing clump 
masses to Jeans masses in TableQ] 

It was difficult to devise an automatic method for identifying 
candidate ETs in the 3D simulations. Instead we generated images 
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Figure 2. Volume Rendering of a 3D simulation of the photo-ionisation of a random clump distribution. Neutral gas density is plotted on a log scale as 
indicated (in units of cm -3 ), with gas at njj < 10 4 cm -3 made transparent to show the difference between neutral tails and dense pillars. Figures are shown 
for outputs every 100 kyr from to 500 kyr. The source is 2 pc off the top of the domain. Note the pillar labelled 'A' front and centre in the third panel, and 
another labelled 'B' at right back in the later panels. 



of the volume renderings from a number of viewing angles and 
picked out structures with n ^ 10 4 cm -3 , an aspect ratio > 2, and 
which persisted for > 50 kyr. Any candidates were then studied 
more closely with cross-section and isosurface plots to verify that 
they could be classified as ETs. 

3.3 Large 3D Simulation 

Based on these tests we tried an intermediate model (Model 16 
in Table fT) with less extreme initial densities. We chose n c i — 
4 x 10 3 cm -3 , intermediate between two models in the tests, with 
clumps restricted to the central 75 per cent of the domain in y and 
z directions, and to 45 per cent of the longer x domain. This gave 
a mass in clumps of 266 Mg, making up ~ 90 per cent of the total 
gas mass. This was distributed among 40 clumps ranging in mass 
from 2.4 — 10.4 Mq, and in radius from 0.06 — 0.14pc. The peak 
density in the initial conditions was tih — 1.8 x 10 5 cm -3 , and 
the background density was rib = 100 cm~ 3 . The ionising source 
was placed 2pc off the domain in the long direction and had an 
ionising flux of 10 11 cm~ 2 s~ x at the front edge of the grid. This 
model successfully produced three ET-like structures and we re- 
peated it with 512 x 384 2 grid cells with a similar physical size 
of 2.67 x 2 x 2 pc 3 . A volume rendering of the neutral gas num- 
ber density is shown in Fig. [2] at 100 kyr intervals from the initial 
state to 500 kyr. The number density is on a log scale, with gas at 



n H < 10 4 cm -3 made transparent (the transfer function is shown 
beside the log scale). This highlights only the densest structures to 
distinguish between shadowed regions and dense tails. 

After 100 kyr the front clumps have begun to accelerate away 
from the source and have clearly been compressed significantly. 
The first pillar (labelled 'A') forms near the front of the domain af- 
ter about 170 kyr and lasts for about 50-70kyr before being acceler- 
ated back to the main body of neutral gas. It is about 0.5 pc in length 
with an aspect ratio of 2 — 3. This frame also shows a large mass 
of merged clumps at the far right. By the next frame at 300 kyr, 
this mass has begun to look more like a pillar, along with another 
mass at the far left. At 400 kyr both of these structures clearly re- 
semble ETs, with the one on the right (labelled 'B') slightly more 
developed and measuring about 0.75 pc in length. All neutral gas 
coloured yellow to red has a number density nu > 10 4 cm -3 , so 
these are dense enough to be considered pillars and not just shad- 
ows. The pillars are formed almost entirely from dense clump ma- 
terial that has been pushed into the shadowed regions. In the last 
frame at 500 kyr pillar B has a surprisingly similar morphology to 
the largest pillar of the three in Ml 6, leaning over its neighbours. 
This pillar has survived for 200 kyr in the simulation, and is likely 
to live for significantly longer. 

We tested the numerical convergence of these results by com- 
paring the structures that formed when this calculation was per- 
formed at resolutions of 128 3 and 256 3 (with shorter x domains) 
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Figure 3. Volume Rendering of a 3D simulation of isolated pillar formation with three massive clumps initially almost in a line, using cooling model CI. 
Neutral gas number density is plotted (cm -3 ), with low density gas made transparent to show the difference between low density neutral tails and dense 
pillars. The transfer function is shown beside the log scale to the left. Simulation times are indicated on each panel. 



with the simulation above. Positions, sizes, and lifetimes of the pil- 
lars were found to be consistent at all resolutions. Transient peak 
densities increased with resolution, as expected in problems involv- 
ing strong compression. 

We summarise the basic results from these 3D models as fol- 
lows. It is easier and more common for a configuration of clumps 
to produce a structure like a mountain or corrugation than an ET. 
In particular, in our 3D models we required clumps of a few solar 
masses in a density field that had substantial low-density regions 
surrounding the clumps in order to generate ET-like structures. Our 
best model produced a single short-lived clump after about 170 kyr, 
and two or three longer-lived pillars after 300 kyr. No models pro- 
duced an ET in less than 150 kyr. The pillars we formed were 
0.5 — 0.75 pc long, had number densities of > 10 5 cm -3 at the 
densest part of the head, and of 1 — 4 x 10 4 cm -3 in the denser 
parts of the trunks. These properties are similar to observed ETs. At 
400 kyr, there was 110 Mq of neutral gas in the domain, of which 
65 Mq had nu > 10 4 cm -3 , divided among two ETs and the other 
clumps. This gives somewhat lower masses in our ETs than in the 
M16 pillars (masses ~ 30 — 90 Mp), but com parable to estimates 
for other observed ETs (e.g. lGahm et alj|2006h . 



4 ANALYSIS OF CLUMP CONFIGURATIONS 

Certain configurations of clumps which occurred in the random 
initial conditions gave rise to pillar-like structures. To investigate 
how these structures are formed, we have simulated configura- 
tions of clumps in isolation to model the formation of pillars A 
and B. These are models 17 and 18 in Table [T] The calculations 
used a computational domain with 192 x 128 2 zones covering 
2.25 x 1.5 2 pc 3 , a source at — 2pc in the x-direction with a lu- 
minosity in ionising photons of L 7 = 2 x 10 50 s _1 and an ambient 
gas density of 200 cm -3 . We have also run these simulations with 
49 small random clumps superimposed on this, corresponding to 
an extra mean density of 400 cm~ 3 , but found very little differ- 
ence in the results. Superimposed on this density field are three 
massive clumps of up to 28.3 Mq each, with a peak density of 
up to Tin — 10 s cm -3 (overdensity of 500) and Gaussian pro- 
files with scale radius 0.09 pc. We have scaled up the mass and 
radius of the clumps compared to the previous simulations in order 
to get closer to the masses and sizes of the M16 pillars, which are 



~ 30 — 90 M© . As expected for simulations without self-gravity, 
the results with less massive clumps were very similar except that 
the time-scales were shorter. In the results that follow we vary the 
relative positions of the three clumps. The ionising flux at the front 
side of the clumps is about 2.9 x 10 11 cm" 2 s -1 . The gas is ini- 
tially neutral, with constant pressure Pg/ks — 2 x 10 s cm -3 K, 
corresponding to 1000 K in the lowest density gas. 



4.1 Three almost co-linear clumps (Pillar B) 

We first investigate the effect of several clumps partially shad- 
owing each other in a roughly linear fashion (model 17), simi- 
lar to the configuration that formed pillar B. The three massive 
clumps are located in the same y-plane as the source at positions 
[0.30,0.72,0.72], [0.75,0.72,0.84], and [1.2,0.72,0.60] (mea- 
sured in parsecs from a corner of the domain nearest the source). 
The front clump thus partially shadows the two behind it. 

Volume renderings for the run using cooling model CI are 
shown in Fig. [3] at times 0, 50, 150, and 250 kyr. At later times 
the clumps merge and the structure no longer resembles a pillar. 
Projections at angles of 20° and 40° to the z-axis (where 0° is 
perpendicular to the pillar) are shown in Fig. [4] with column den- 
sity shown in the maps, and line of sight (LOS) velocity shown in 
position- velocity (P — V) diagrams below the maps. The 3D P — V 
datacube is projected on to 2D by summing the contribution of all 
the pixels in a given image y-column at each x position. This means 
the normalisation of the P — V column densities is on a somewhat 
arbitrary scale. 

Fig.[3]shows how the partially shadowed clumps compress and 
move sideways into the shadow, while the fully exposed clump is 
accelerated and merges into them. In this situation the gas is ini- 
tially aligned as a pillar-like structure, and the compression due to 
photo-ionisation of the surroundings serves to enhance this appear- 
ance for a limited time. The neutral gas only resembles a pillar for 
150 kyr; this structure is a relic from the initial conditions as op- 
posed to being generated by dynamical evolution. The P — V di- 
agrams in Fig. [4] clearly show gas moving at very different speeds 
away from the radiation source, so it is likely that some of this gas 
will move into the tail as the structure recedes from the source and 
is subjected to lower ionising fluxes. 

Even though model 16 did not show such a clear early-forming 
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Figure 4. Projections through the "co-linear clumps" 3D simulation shown in Fig.[3]with cooling model CI. Gas column density of neutral (atomic) hydrogen 
is plotted in a 2D projection above and a position- velocity diagram below (with radial velocities indicated in units of kms" 1 ). Times are indicated on each 
panel. The normal vector for the projections is at angles of 20° and 40° from the z-axis in the left and right sequences respectively, and is perpendicular to the 
y-axis. Column density is shown on a log scale as indicated (in units of cm -2 ). 



pillar, it is nevertheless interesting that the projected column den- 
sities in Fig. |4] show structures that resemble pillars for the first 
150 kyr of the simulation. The clumps were not as neatly arranged 
in the random clumps model, but given that the ISM in molecular 
clouds tends to be filamentary and clumpy, it is certainly possible 



that there are some ETs in H II regions which are formed in this 
way i.e. purely from initial conditions and not dynamically. This is 
a less satisfactory explanation for Ml 6 however, since we would 
need three lines of clumps/filaments positioned beside each other, 
and all pointing by chance back to the brightest star in the nebula. 
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Figure 5. Volume Rendering of a 3D simulation of isolated pillar formation with three massive clumps initially in a triangle configuration, using cooling model 
CI. Neutral gas number density is plotted ( cm -3 ), with low density gas made transparent to show the difference between low density neutral tails and dense 
pillars. 



Looking at the evolution in more detail, the P ~ V diagram 
shows the acceleration of the first clump away from the source, up 
to about 6 kms -1 in 150 kyr (The true velocity is v\\/ s'm6). For 
all of this time the shadowed gas remains essentially stationary. The 
velocity signature of this formation mechanism is very clear, with 
the head of the pillar receding from the star faster than the trunk. 
This is the opposite o f what is seen in o t her fo rmation scenarios 
such as the model of iLefloch & LazarefJ J 19941) of a tail forming 
behind a single clump which has the tail streaming away from the 
star faster than the head. We show the two perspectives at 20° and 
40° to demonstrate that even though 20° is close to perpendicu- 
lar, they both show the same trends in velocity along the pillar (In 
subsequent figures we only show 20° projections). 

We note, however, that the almost linear decrease in LOS ve- 
locity from the head to the tail is significantly enhanced by trans- 
verse motions in the gas. The middle clump is being pushed away 
from the observer into the shadow, and the right-most clump to- 
wards us (cf. Fig. [3](- If we make the same projections from the 
opposite side of the simulation these transverse velocities are re- 
versed leading to a much less obvious velocity gradient (although 
the large velocity at the head is still evident). We return to this issue 
later in Section [5] where we show velocity profiles from different 
perspectives. 



4.2 Triangle of three clumps (Pillar A) 

We now investigate a more dynamical situation (model 18), in 
which two clumps are pushed past a third which was initially shad- 
owed, and are wrapped around into the tail, forming a pillar. This 
was how the short-lived 'pillar A' formed in model 16, and it also 
crudely models the later evolution of pillar B. In this scenario, 
the three massive clumps are placed at positions [0.45, 0.72, 0.6], 
[0.45,0.72,0.9], and [0.75,0.72,0.72], giving us two clumps in 
front shadowing the third clump, all in the plane y — 0.72 pc. This 
plane also contains the source at [—2.0,0.72,0.72], The proper- 
ties of the two front clumps were modified to have masses M — 
27.7 Mq, peak overdensities of 250, and scale radii of 0.1125 pc, 
while the shadowed clump has mass of M = 28.3 Mq, overdensity 
of 500, and scale radius of 0.09 pc, as in the previous model. The 
formation of a pillar in this configuration is quite sensitive to the 
relative clump positions; if the front clumps are too close together 



or too dense they do not move past the shadowed one without dis- 
rupting it and if they are too far apart they never merge with the 
shadowed clump to form a single clump/pillar configuration. 

We first ran this simulation using cooling model CI and a 
volume rendering of the results is shown in Fig. [5] The pillar de- 
velops more slowly than in the previous scenario, taking around 
200 — 250 kyr to form, because the gas has to travel further to 
get into the tail. Initially the front two clumps are compressed 
and slowly accelerated, seen after 50 kyr in Fig. [5] They collide 
obliquely with the shadowed third clump, which then starts to ac- 
celerate due to a combination of being exposed to ionising radiation 
and to the shocks driven through it by the passage of the first two 
clumps. Material from the remains of the first two clumps sweeps 
into the tail region after 200 kyr, producing the pillar-like structure 
seen in the last panel at 250 kyr. In this model the ET structure is 
formed dynamically rather than being left over from initial condi- 
tions, and the column of neutral gas resembles a pillar only after 
250 kyr, unlike the previous model. 

As one might expect from the initial conditions, this column is 
quite asymmetric; it is much broader in one direction than the other 
at late times. The volume rendering shows it from the narrower 
perspective; projections at 20° from perpendicular to the pillar are 
shown in Fig.|6]with edge-on projections at right and face-on at left. 
The gas in the trunk is at volume densities of 1 — 3 x 10 4 cm -3 , 
but projecting through the narrow axis does not build up enough 
column density to resemble a pillar whereas the projection through 
the broad axis clearly resembles an ET at 250 kyr. 

The velocity profiles show a broad range of velocities which 
vary significantly with perspective and over the evolution of the 
system. We do not see a strong signature here as was seen in Fig. [4] 
One interesting feature is seen in the right P — V diagram at 
150 kyr. This shows high velocity gas from the front two clumps 
which has pushed past the third and moved i nto the trunk away 
from t he head. This is a similar profile to that of lLefloch & Lazarefj 
( 1994) for gas streaming from the head into the tail of a cometary 
globule. The feature is not clearly seen in the left-hand panel be- 
cause the tail is strongly re-expanding along the LOS, leading to 
a broad velocity profile which has the opposite slope from head to 
tail compared to the right-hand panel. Here transverse motions are 
once again masking the true trend in recession velocity from the 
star along the pillar. 
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Figure 6. Projections through the "triangle of clumps" 3D simulation shown in Fig.|5]with the CI cooling model. The projection is at 20° to the y-axis on the 
left, and to the 2-axis on the right. The two front clumps are projected on to each other in the right sequence of figures. 



5 EFFECT OF COOLING ON PILLAR FORMATION 



iMellema et all {2006a) use a cooling model similar to our CI 
model, with relatively little neutral gas cooling, to model the global 
expansion of an H II region into a turbulent density field. Individ- 
ual features are not highly resolved, but it is clear that their results 
(e.g. their fig. 5) are qualitatively similar to ours, showing fairly 



smooth coherent pillar-like structures. In particular both they and 
we do not find I-front instabilities developing, or any fragmenta- 
tion which is a feature of simulations with strong neutral cooling. 
This is u ndoubtedly due t o the thermal physics, since it has been 
shown bv lWilliamsi d2002h and lWhalen & Normanl d2008t) . among 
others, that instabilities in D-type I-fronts are strongest when the 
neutral gas can cool rapidly. These instabilities can generate dense 
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Figure 7. As Fig. [3] but with the simulation run with cooling model C2. Simulation times are indicated on each panel. 
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Figure 8. As Fig. [3] but with the simulation run with cooling model C2, and with frames at different times. 



finger -like structures (e.g. IWhalen & N orman 2008; Mizu ta et al.l 
2006) but it is unclear if they could generate something as massive 
and large as the Ml 6 complex of pillars. 

Other authors (e. g. iLefloch & Lazarefdll994 1 Willia ms et al.l 
l200ll : lLora eTai]|2009l ; lGritschneder et al.ll2009h have used a ther- 
mal model where the neutral gas is cold and isothermal, leading 
to st rong fragmentation and instability in the I-front. Observation- 
ally, IWhite et alj dl999h showed that the Ml 6 pillars have cold 
cores (T ~ 10 — 20 K), warm trunks or fingers (T ~ 60 K), and 
a 'hot' shell of surrounding gas at T ~ 250 — 3 20 K. To study the 
possible effects of these temperature variations, iMiao et al.l (2006) 
use a much more detailed thermal model to study the head of a 
pillar numerically, but at the expe nse of using very low resolution 
for the dynamics. More recently, iHennev et alj d2009l) have fit a 
number of functions to detailed calculations of cooling and heat- 
ing rates for their MHD models of photo-evaporating clumps. They 
found the shadowed neutral gas was not isothermal, but had a range 
of temperatures comparable to observed values. These results show 
that while the isothermal approximation is expected to be more re- 
alistic than our CI model, it is unclear whether it is sufficient to 
capture the details of ET formation. 

To investigate the influence neutral gas cooling has on our re- 
sults, we have repeated the simulations in the previous section using 
cooling model C2, designed to be intermedia te between CI and th e 
isothermal approximation. In comparison to lHennev et alj j2009t) . 



their 'molecular' cooling rate (which they denote Lpdr) scales as 
L oc p 16 erg cm -3 s _1 whereas our C2 model scales as L oc p. So 
while our model has significant neutral gas cooling down to 100 K, 
very dense gas cools more rapidly in their simulations, closer to an 
isothermal model. 

Volume rendering of the linear and triangular clump configu- 
rations with C2 cooling are shown in Figs.|7]and[8] The evolution 
is slower than with CI so we show the results at times 50, 150, 
250, and 350 kyr on the same neutral gas number density scale. 
Comparing the 50 kyr panels between the CI and C2 runs, the 
main change is immediately apparent. The extra neutral gas cooling 
makes shocks more compressive and slower moving. This can be 
clearly seen in the shadowed region, where the converging shock is 
about 10 x denser in Fig.|7]than in Fig. [3] The cooling also affects 
the exposed clump; the radiatively driven implosion phase is slower 
and the clump is compressed to a greater degree. The rocket effect 
is less effective on this exposed clump: firstly there is a smaller sur- 
face area to absorb the photon flux; secondly, the denser gas in the 
I-front and in the photo-evaporation flow leads to faster recombi- 
nations and thus more photons are required to ionise an atom and 
keep it ionised until it flows away from the I-front. 

Comparing the 150 kyr and 250 kyr panels of Figs. [3] and [7] 
the difference is dramatic. The leading clump has fully merged into 
the second in the CI run, and the implosion phase is followed by 
a strong re-expansion. Re-expansion in the C2 run is much weaker 
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Figure 9. Projections through the "co-linear clumps" 3D simulation shown in Fig.[7]wifh cooling model C2. Only the projection at 20° to the z-axis is shown, 
but from the front of the domain at left and from the rear at right. This highlights the effect transverse gas motions can have on velocity profiles. 



because most of the heat generated is radiated away before it can 
instigate a rebound. We obtain a much narrower and denser struc- 
ture which resembles an ET for at least 200 kyr. The initial config- 
uration is enhanced rather than disrupted by the photo-ionisation 
process. 

Fig. [8] shows a very similar story for the triangle of clumps. 
The system's evolution is slower, the shocks are more compressive, 



and there is less re-expansion after the implosion phase. The gen- 
eral picture from these two models is that ETs are longer, narrower, 
denser, longer-lived, and have more substructure with the C2 cool- 
ing prescription. While qualitatively the same evolutionary scenario 
plays out, quantitatively the neutral gas cooling has a significant ef- 
fect on the results. 

Projections through the C2 runs are shown in Figs.[9landllOI 
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Figure 10. Projections through the "triangle of clumps" 3D simulation shown in Fig.[8]with the C2 cooling model. The projections are at 20° to the y-, 
the left, and to the z-axis on the right, as in Fig.|6]but shown at different times to reflect the slower evolution of this model. 



Fig. [9] shows the projection at 20° as in the left sequence of Fig. [4] 
for CI, but the two sequences show projections from the front and 
rear perspectives. The column density maps are similar from both 
perspectives and it is clear that the structure resembles an ET at 
all times shown. The left P — V diagrams show the same trend as 
in the CI model, where gas at the head is receding more rapidly 
from the star than the trunk is, but this trend takes significantly 



longer to become established. The right hand P — V diagrams also 
show the head receding rapidly, but the gas from the partially shad- 
owed clumps has very different LOS velocities due to transverse 
motions. This particular case shows how difficult it is to infer any- 
thing about a formation scenario from aP — V diagram. With non- 
axisymmetric initial conditions, transverse motions mix with the 
recession velocity in an unpredictable way and are a strong con- 
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taminant when one wants to measure the recession velocity along 
the length of a pillar. Fig.llOlshows projections face-on and edge-on 
through the triangle of clumps scenario, again at 20° to perpendic- 
ular. The same asymmetry as shown in Fig. [6] is apparent, but to a 
lesser degree, and the C2 run shows that the column of gas resem- 
bles an ET in both directions from 200 — 300 kyr, albeit much more 
convincingly in the left hand sequence. The P — V diagrams show 
very little trend in velocity along the length of the pillar at most 
times, but again it can be seen in the 200 kyr panels that there is 
high velocity gas moving from the head into the tail region, as was 
seen in the CI model. 



6 DISCUSSION 

Willi ams et ail feOOlb performed axisymmetric simulations with a 
similar aim to our work - to investigate the mechanisms by which 
ETs can form. They modelled parallel rays rather than treating the 
rad iation from a source a t a finite distance. In most of their mod- 
els, l^lliamseraL I feOOlh start with a dense (n ~ 10 cm ) layer 
of gas at the boundary most distant from the star. This served both 
to stall the I-front and to provide a reservoir of dense gas from 
which to build up a pillar. Th is is a very similar picture to that mod- 
elled more recently in 3D by Ra ga et al] d2009l) . who studied small 
(I ~ 0.03 pc) columns of dense gas forming in a photo-ionised re- 
gion. Their initial conditions had an effectively infinite reservoir of 
neutral gas shadowed by smaller clumps. Evaporating gas from the 
reservoir flowed into the shadowed regions forming dense columns 
as it recombined and cooled. These authors have demonstrated a 
mechanism by which ETs can form given a sufficiently large reser- 
voir of dense gas behind a shadowing clump. We consider alternate 
scenarios, however, where the ET must be built up from clumps of 
comparable mass and from low density inter-clump gas. We find 
that this is sufficient to produce pillars when the gas is in certain 
configur ations, but we also find that they do not form as readily as 
found by Willia ms et alj j200 lh . This is quite important: while ETs 
are seen in many H II regions, structures resembling heaps and 
corrugations which are not elongated are much more common. The 
mechanisms by which ETs form cannot therefore be too efficient 
or many more ETs would be observed. 

The fact that most H II r egions seem to have one or two ET 
structures led IWilliams et alj J200 lh to suggest that they may be 
long-lived objects, with lifetimes comparable to the H II region, 
a proposition supported by their simulation results. An alternative 
scenario is that they are short lived (~ 100 kyr), but multiple gen- 
erations of them occur during the expansion of an H II region. For 
models where the gas is not already organized in a linear structure, 
we do not find any pillars forming in less than 150 kyr. Adding extra 
cooling in the neutral gas only lengthens this formation time. Addi- 
tionally, it takes about 250 — 300 kyr before more massive and long- 
lived pillars start to form dynamically. This is the rough time-scale 
for dense neutral gas to be accelerated, pushed past other dense 
clumps, and stretched out into a long tail. We suggest it would 
be difficult to dynamically generate a parsec-scale dense ET from 
static initial conditions in less than 150 kyr, and our results are con- 
sistent with the claim bv lWilliams et alj ( 12001!) that ETs are likely 
to be long-lived structures. 

6.1 Gas Pressure 

Observational age constraints are thus far not very stringent. An 
upper age limit for ETs is set by the age of their H II region. 




Figure 11. Simulated H a emission map from our simulation of a pillar 
forming from three almost co-linear clumps using cooling model C2. The 
map can be compared to projected density in the bottom left panel of Fig. [9] 
The scale used is linear with the strongest emission in black and deliberately 
saturated to bring out low contrast features. 

The free-fall time of the densest parts of the pillars does not con- 
strain their age as the gas is being actively compressed and so 
the free-fall time is changing as the structures evolve. If the pres- 
sure difference between the neutral and ionised gas is very large, 
however, this may indicate that the dense gas has not had time 
to dynamically respond to the photo-ionisation, and so the struc- 
tures must be young, perhaps younger than their sound crossing 
time. In M16, the gas pressure at the b ase of the ionised photo- 
evaporation flow has been estimated by iHester et alj d 19961) to be 
Pg/ks — 6 x 10 7 cm -3 K. The inter nal pressure in the dens- 
est parts of the pillars was estimated bv lWhiteetal]dl999h to be 
Pg/ks = 3.5 x 10 7 cm -3 K. These values are based on temper- 
ature and density e stimates in ionised and neutral gas respectively. 
IWhite et alj J 19991) interpreted this pressure difference as suggest- 
ing the ETs in Ml 6 may be young, with the dense clumps currently 
under going the early stages of radiation-driven implosion dBertoldil 
1989). In our simulations the peak pressure at the ionisation front 
varies between Pg/ks = 5 x 10 7 and 2 x 10 8 cm -3 K, and in 
the neutral ETs between Pg/ks — 2 x 10 6 and 5 x 10 7 cm -3 K. 
These values are consistent with the observations, but we found that 
an equilibrium state was never reached in our models. The pressure 
varies by large factors both in time and in space along the length 
of the ET, and the peak pressure changes in time according to the 
instantaneous density at the I-front. These results suggest that the 
observed pressure difference (a factor of 1.7) should not be taken 
as strong evidence that the pillars are young since we find pres- 
sure variations much larger than this within the dense neutral gas. 
A caveat to this is that our thermal modelling of the neutral gas is 
crude and gas pressure is sensitive to this, so it would require more 
detailed modelling to make a definitive statement. We can say that 
simulations with both cooling models CI and C2 have these pres- 
sure variations, and we believe the dynamical nature of the ETs 
which form in our models will always generate significant pres- 
sure gradients within the ETs. It would be interesting to compare 
this with structures formed via I-front instabilities to see how their 
internal dynamics differ. 

6.2 Morphology of ETs 

Molecular emission traces the projected mass density of ETs and 
has show n them to be rather more clumpy than they appear in op- 
tical data jPoundl 19981 ; IWhite et al.ll 19S*gf) . with significant density 
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variations along their length. The heads of the pillars are the dens- 
est regions, with clumpy lower density gas found in the trunk. Our 
simulations support this clumpy model for the underlying density 
structure; in some of our models the clumps maintain their iden- 
tity for hundreds of kyr (e.g. Fig. [9j. We also find that the high- 
est density gas is always found in the head of the pillars because 
the strongest compression is always found ahead of the I-front 
which drives the dynamics. We have also calculated optical emis- 
sion maps due to recombination ra diation (e.g. H a ) in the man- 
ner d escibed by Henn ev"et alj ^20051) . using a constant dust opacity 
as in Melle ma et alj j2006a ) . An image from the co-linear clumps 
simulation (model 17) at 250 kyr is shown in Fig. [TTJ this can be 
compared to the bottom left panel of Fig. [9] which shows projected 
neutral gas density and is closer to what we expect from a molecu- 
lar emission map. It is clear that the 'optical' image shows a much 
smoother structure which appears more like the HST images of 
M16. In particular the dense parts of the pillar's trunk are com- 
pletely dark in the optical image while there is substructure in the 
projected density. Our image has less substructure than the HST 
images; this is likely due to the limited spatial resolution of our 
simulations, and possibly because our cooling models do not pro- 
mote the development of I-front instabilities. 



6.3 Position-Velocity Diagrams 

The full 3D nature of our simulations has allowed us to calcu- 
late line-of-sight velocity profiles through the ETs which include 
asymmetric transverse motions. We have shown that these motions 
can largely determine the observed profile (Fig. [9}. This has im - 
plic ations for in t erpreti ng the P — V diagrams in IPound Jl998t) 
and IWhite et alj <1999b . in particular for the largest M16 pillar 
which has a velocity gradient that changes sign from the head to 
the tail. This was interpreted as evidence that the head and tail are 
separate structures, but this is difficult to reconcile with optical ob- 
servations. Our results suggest a resolution to this apparent con- 
tradiction: this pillar could have started out as two separate dense 
clumps which have since merged but have kept their identity in 
velocity-space. Fig.[9]shows just such a case where the clumps gain 
opposing transverse velocities, giving the projected velocity profile 
an S-shape which is maintained for > 100 kyr despite the clumps 

having merged into a pillar much earlier. 

For the smallest of the M16 pillars. IPound j 19981) notes that 
the shadowing implies it is closer to us than the ionising stars so we 
expect the tail to stream towards us faster than the head. In fact the 
opposite is true. Our results also offer two possible explanations 
for these observations: (A) it could be due to a strong transverse 
motion in the pillar if it is seen close to edge-on, or (B) it could be 
formed from a column of dense gas closely aligned with the radi- 
aton propagation direction. Our models of this situation show that 
for at least some of the pillar's evolution the head can be receding 
from the star faster than the tail. 



6.4 Influence of other physical processes 

We have already discussed in detail the effects of neutral gas 
cooling on our results, but there are a number of other pro- 
cesses which could co ntribute significantly. Since mass ive stars 
generate strong winds, iRaga. Steffen. & Gonzalez! 120051) studied 
photo-evaporating clumps interacting with a stellar wind, finding 
that when the photo-ionisation was sufficiently intense the photo- 
evaporation flow effectively shielded the clumps from the wind. 



This gives us co nfidence that ion ising radiation is the main driver 
of ET formation. IRaga et alj (2009) studied photo-ionisation mod- 
els with a basic treatment of diffuse radiation. They found there 
was no strong difference between runs with the diffuse radiation 
and runs where they used the on-the-spot approximation. This sug- 
gests that, while diffuse radiation would undoubtedly change our 
results somewhat , it i s unlikely to make a dramatic difference. 
lEsauivel & R aga (2007) studied the effects of self-gravity on the 
photo-ionisation of a dense cloud of gas using a two-temperature 
equation of state (i.e. the neutral gas is isothermal) which gives 
rise to strong instabilities and fragmentation associated with the 
photo-ionisation process. Interestingly, they found that self-gravity 
had little effect on the overall process of the evaporation and frag- 
mentation of the massive clump, and was only significant in de- 
termining the properties of the densest sub-clumps produced by the 
fragmentation. This supports our implicit assertion in this work that 
pressure forces are the dominant driver in the evolution of photo- 
evaporating clumps, at least in the early stages. 



7 CONCLUSIONS 

Using moderately high resolution 3D radiation-hydrodynamics 
simulations of clumpy density fields exposed to ionising radiation 
from a point source, we have investigated how effectively shadow- 
ing can generate pillar-like structures. We have shown that even in 
the absence of self-gravity or I-front instabilities, large parsec-scale 
ETs can form dynamically solely due to this shadowing. 

Of our 2D and 3D models, about 20 per cent produced long- 
lived pillars. This is more due to our choices of initial conditions 
than how easy or difficult it is to generate ETs. Nevertheless we did 
find certain density fields more conducive than others to forming 
pillars. The most successful of these contain both large low density 
regions where the I-front can propagate far from the source, and 
also massive clumps with sufficiently high density and inertia to 
stall the I-front for > 100 kyr. Pillars are formed with diameters 
comparable to those of the clumps that give rise to them. 

Simulations of specific configurations of massive clumps 
show that pillar-like initial conditions evolve with a different ve- 
locity signature to configurations where the ETs form dynamically 
from clumps that are not initially co-linear. For a simulation with 
three clumps initially almost co-linear we find the head of the ET 
recedes from the source more rapidly than the shadowed trunk 
which has not been exposed to radiation. For dynamically form- 
ing ETs, gas streams into the shadowed trunk past the head, and is 
thus moving faster than, or at a comparable speed to, the pillar's 
head. This could produce a noticeable observational signature from 
which these two formation mechanisms can be distinguished, but 
we find that the P — V diagrams vary significantly with viewing 
angle due to transverse gas motions. This variation increases with 
the degree of asymmetry and the transverse motions significantly 
contaminate attempts to measure recession velocity in the pillars. 
These transverse motions offer an explanation for the features seen 
in P — V diagrams for the pillars in Ml 6. 

We have shown that neutral gas cooling has a very strong in- 
fluence on the modelling results. Stronger cooling produces ETs 
which take longer to form dynamically, are narrower and denser, 
and are more resistant to the rocket effect and hence live longer. 
This shows that the complex chemistry and thermal physics in 
molecular clouds may play a crucial role in ET formation and evo- 
lution. For the specific case of almost co-linear clumps, the initial 
pillar-like configuration was slowly disrupted in simulations with 
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the CI cooling model. By contrast it was enhanced with the C2 
cooling model, also proving to be long-lived. 

We have not addressed magnetic fields in this work. 
lHennevetallj2009h have studied the photo-ionisation of a magne- 
tised globule, finding that strong uniform fields can have a signifi- 
cant influence on the evolution of the photo-ionisation process. In 
future work we will investigate the effects of more realistic thermal 
physics, as well as the presence of magnetic fields, on the results 
we have presented here. 
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